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ABSTRACT 



A digital computer program which will provide accurate and stable 
solutions of linear and nonlinear ordinary differential equations is 
described, as are subprograms which simulate many of the nonlinearities 
occurring in feedback control systems. The program employs fourth 
order Runge-Kutta numerical integration with automatic error checking 
and interval modification. Specialized coding sheets assist in the 
preparation of input data, and built-in print and graph output routines 
provide a permanent record of input equations^ parameter values and 
output data for each solution. The program is flexible, yet it can 
be used by a person with only a rudimentary knowledge of FORTRAN 
programming. 
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1 . 



Introduction. 



The analysis of scientific problems frequently leads to a set of 
simultaneous ordinary differential equations which describe the phenomenon 
being studied. Analytic solution of these equations is often impracti- 
cable and analog or digital simulation becomes necessary. A digital com- 
puter program which will solve ordinary differential equations has been 
written in FORTRAN 60 for the CDC 1604 Digital Computer. It includes 
FUNCTION subprograms which cover many of the nonlinearities occurring 
in feedback control systems [^7] . The program has been named SIDES j, 
an acronym for Simultaneous Differential Equation ^olver. 

Digital simulation should be considered as a valuable supplement 
for the analog computer, but not as a replacement. Selection between 
analog and digital simulation methods should be based on an evaluation 
of the problem to be solved: the complexity of the equations^ the magni- 

tude of the included parameters and the required accuracies of the solu- 
tions. Factors which favor the selection of digital simulation^ and 
more particularly, SIDES, are: 

a. The necessity for number scaling is eliminated. 

b. Solution accuracies are achieved which are unattainable 
with an analog computer. 

c. Printed and graphical outputs provide a permanent written 
record of input equations, initial parameter values and output data for 
each solution. 

d. Nonlinear operations such as squaring^ forming trigono- 
metric and exponential functions, and multiplication are easily performed. 
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e. Special functions such as hysteresis and time delay are 
available or can be readily generated as ordinary FORTRAN SUBROUTINE 
or FUNCTION subprograms. 

Some of the limitations ares 

a. Inherent errors are due to the sequential operation of the 
digital computer in contrast to the simultaneous processing of signals 
perfoxmied by the analog computer. 

b. The effect of parameter variations on problem results cannot 
be immediately observed. (This limitation could be partially eliminated 
by the incorporation of a data display unit with the CDC 1604). 

Existing digital simulation programs vary considerably with regard 
to the input language and the integration techniques used. Stein^ Rose 
and Parker [ 5 ] have devised a system which accepts as input the descrip- 
tion of a ’’patching** diagram prepared for the Electronic Associates PACE 
Analog Computer. FORTRAN is used as an intermediate in the compilation. 

The DYSAC program written by Hurley [2] is in symbolic machine language 
(CODAP). It contains a supply of digitally simulated analog computer 
components which are interconnected by the ’’patching” information sup- 
plied by the user. Integration is performed by the fourth order Runge- 
Kutta method with a fixed Integration interval. Nordsieck [4] has 
developed a method for automatic digital computer integration which is 
equivalent to a reformulation of the method of Adams. 

The SIDES program to be described here^ was designed with particular 
emphasis on user convenience and solution accuracy. Input and output 
formats are built-in and integration is performed by the fourth order 
Runge-Kutta method with automatic error checking and interval modification 
as developed by Anderson [l^ . A listing of the SUBROUTINE and the FUNCTION 
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subprograms which have been developed for use in conjunction with SIDES 
is given in Appendix B. 
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2. The SIDES Input Routine. 

Calling program . SUBROUTINE SIDES will normally be called from the 
FORTRAN Library tape by a calling program which must be of a prescribed 
form. To simplify the writing of this program, a SIDES Coding Sheet has 
been prepared (Appendix A) which reduces most of the job of programming 
to the simple task of filling in blanks and lining out undesired options. 
The calling program deck must be of the following form; 

..JOB* JONES J BOX 113 (Job ident.) 

PROGRAM PROBLEM (Arbitrary name) 

DIMENSION X(50), XDOT(30), C(50) (Standard card) 

C(50) = 1.0 

1 CALL SIDES (T,X,XD0T,C) 

A 

V 

GO TO 1 
END 
END 

This deck is followed by a blank card and the data and option cards. 

A typical program, together with its output, is included in 
Appendix A. 

Equation statements . While equations for different problems vary, 
they can be manipulated into the standard first order form; 



ft ff 

t9 99 

(Equations written 
in FORTRAN language) 
(Standard card) 

99 19 

19 II 



i = f. (x, jX,, 

dt 1 i ^ 



X ,t) 

n 



r 1,2,. ..,n 



in which the independent variable must be designated by t and the 
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dependent variables by The number of such equations is designated 

by n with n 30. Since only the x and t variables can be output ^ addi- 
tional x^, n < i^ 48, have been made available. If some quantity other 
than one of the variables is desired as output, it must be equated to 
one of these x^. 

A set of constants, c^, 1"^ i^ 49, has also been made available so 
that parameters may be introduced or varied by data cards. 

Any legal FORTRAN statement or technique may be used in writing 
these equations, provided that a transfer of control does not prevent the 
processing of all relevant statements during each step of numerical 
integration. All of the FORTRAN Library functions are available for 
appropriate inclusion in the equation statements. Some additional 
FUNCTION subprograms which have been developed for use by SIDES are 
described below. These and other subprograms may be appended to the 
calling program in the conventional manner. They may also be added to 
the library tape if appropriate. 

Data and option cards . The seven sections of input data required 
to specify a problem for SIDES are enumerated in Table I. All sections 
are straight- forward and their purpose fairly obvious. 

If several solutions of the same equation are required, it is 
necessary to specify the number of runs and to supply a data deck for 
each run. For each run after the first, only the cards of sections 2 
through 7 (Table I) are required. Information can be retained between 
runs by choosing the appropriate HELD option. If it is desired to sub- 
divide a long solution into two or more successive runs, the final condi- 
tions of one run may be used as initial conditions for the next by making 
use of the FCS HELD command. 
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TABLE I 



DATA REQUIRED FOR A PROBLEM 



DATA SECTION 


CONTROL OPTIONS 


1. Problem description 


1 — >30 first order differential 
equations 
1 — >*9 runs 


2. Coefficients 


Values zeroed 
Values set in 
Values held 


3. Initial conditions 


Values zeroed 
Values set in 
Initial or final condition 
values held 


4, Time 


Initial and final time 
given 

Values held 


5, Integration step 
size 


Step size given 
Step size computed and 
automatically varied 
to maintain a pre- 
scribed precision 
Step size held 


6. Data print-out 


1-^ 8 variables printed 
l—> 99 steps of integration 
between print-outs 
Data held 
No print-out 


7. Graph output 


1 — ^5 single-curve graphs 
1 — ^5 curves on a single 
graph 

1 — ► 99 steps of integration 
between graph points 
Data held 
No graph output 
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3. The SIDES Output Routine. 

All job identification, together with the run number and a record 
of all input data, is automatically output for each run. A maximum of 
eight quantities may be printed-out during integration. An automatical- 
ly scaled graph output is also available, providing a maximum of either 
five single-curve graphs or five curves on a single graph. 

A maximum of 10,000 integration increments, 300 lines of print-out, or 
900 graph points (per curve) can be generated during any one run. The run 
will be terminated if any one of these conditions is exceeded. 
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4. 



Error Indications , 



The SIDES program will be terminated when all runs called for are 
completed or if an error is detected in the calling program or in the 
data and option cards, SUBROUTINE SIDES contains 38 error stops, each 
preceded by a diagnostic print-out. These error stops refer to diffi- 
culties such as incorrectly punched cards, missing cards and improper 
option selection. All error stops and print-outs inherent to the basic 
FORTRAN Compiler also apply. 
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5, The Integration Method. 

The operation of SIDES is centered around the fourth order Runge- 
Kutta method of numerical integration. The choice of integration algo- 
rithm is necessarily a compromise and Runge-Kutta was chosen primarily 
because of the ease both in starting solutions and in changing the 
integration increment during solutions. The use of the method developed 
by Anderson [ 1 \ to compute the truncation error has eliminated one of 
the main shortcomings of Runge-Kutta as used in other programs [6^ . 

On this basis, the integration interval can be automatically varied to 
ensure stability and to keep the error during each step within prescribed 
bounds, without excessive lengthening of the computing time. The allow- 
able truncation error per unit time for each of the dependent variables 
was chosen to be: 



- A 

10 (variable magnitude ) 
i total time 

-9 

The variation of the step size is bounded by arbitrary limits of 10 

-3 

and 10 times the total time. 

In the event that the user wishes to prescribe a fixed step size^ 
the Runge-Kutta integration can be processed in either one or two 
segments, each with a different step size. 
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6. Accuracy Considerations. 

Since the accuracy of the Runge-Kutta method of numerical integra- 
tion is a function of the size of the integration increment, use of a 
digital computer allows solution accuracies that in most cases far ex- 
ceed existing engineering demands. The following two examples are 
given to indicate the effect of integration interval on solution accuracy. 
The interval was held constant for each run. 

Example 1 ^ 



Differential Equation 


y' = -2xy^ 


(y o ~ ^ > X o • 


Explicit Solution 


y = (1+x^-l 




Integration Interval 


y at X = 5 


Error at x = 5 


0.03 


3.8x10"^ 


4xl0‘^^ 


0.05 


3.8x10"^ 


-9 

3x10 


0.10 


3.8xl0'^ 


4x10 ^ 


2 

Example 2 






Differential Equation 




(y X q= 0) 


Explicit Solution 


y = (1+x)^ 




Integration Interval 


y at X = 5 


Error at x = 5 


0.005 


9.4x10^ (at 
X = 2.9) 


1x10 ^ (at x=2 


0.03 


8.1x10^ 


2xl0'^ 


0.10 


8.1x10^ 


4x10^ 



^Milne, W. E. Numerical solution of differential equations. 
John Wiley & Sons., 1953: pp. 20-21. 

^Ibid., pp. 74-75. 
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Although Example 2 is intentionally unfavorable for solution by 
Runge-Kutta integration, it is obvious that high accuracies can be 
achieved, although at the expense of additional computing time. 

One of the main advantages of the computed step size integration 
method is that solution accuracy and computing time are not seriously 
impaired by problem discontinuities. This is an important factor to 
the engineer who must analyze systems containing switching devices or 
step-wise parameter fluctuations. Solution accuracy in SIDES is deter 
mined by the estimated magnitude (power of ten) of the dependent vari- 
ables, X(j), 1 ^ j ^ n. An examination of Example 3 indicates that 
this estimated magnitude can often vary widely before solution ac- 
curacy suffers. 

Example 3 
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s (s+1) 
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where R is of the form: 
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VARIABLE MAGNITUDE 
SETTING 


C AT 
T = 2 


ERROR AT 
T = 2 


C AT 
T = 4 


ERROR AT 
T = 4 


-3 


3.2 


4x1 o' 


- 


- 


-2 


3.2 


5x10-10 


- 2.0 


2 x 10-10 


-1 


3.2 


.9 

5x10 ^ 


- 2.0 


-9 

3x10 


0 


3.2 


IxlO'^ 


- 2.0 


4x10-0 


1 


3.2 


2xl0'0 


- 2.0 


9xl0'^ 


2 


3.2 


7xlQ-0 


- 2.0 


SxlO'O 


3 


3.2 


7x10-0 


- 2.0 


8x10-0 


4 


3.2 


7x10-0 


- 2.0 


9x10-0 
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7. Special FUNCTION Subprograms. 

One of the stated advantages of digital solution was the relative 
ease with which many nonlinearities can be simulated. The develop- 
ment of subprograms to simulate functions such as hysteresis and time 
delay, however, presented special problems in relation to their use with 
Runge-Kutta integration with variable step size. Since FUNCTION simula- 
tion of these nonlinearities requires past information for its proper 
operation, it was necessary to modify the basic SIDES program so that 
there could be a free interplay of information between it and the 
particular FUNCTION subprogram. Utilizing the fact that the COMMON 
statement, as applied to FORTRAN 60, will cause values to be entered 
in known cells, it was possible to establish this needed correspondence 
between the two separately compiled subprograms. 

FUNCTION subprograms have been written to simulate time delay and 
two and three position relays with hysteresis. These subprograms are 
available for use in conjunction with SUBROUTINE SIDES. A brief descrip- 
tion of each FUNCTION is included in Appendix A and a listing of each 
FUNCTION subprogram appears in Appendix B. 
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8. Recommendations and Conclusions. 



SUBROUTINE SIDES has proven to be an extremely easy program to use. 
The SIDES Coding Sheet enables the user to prepare his program with a 
minimum of reference to additional explanatory material. Due to the 
complexity of the input program, it is believed that the coding sheet 
is almost a necessity. Continued use of this sheet will probably point 
out needed changes to increase its effectiveness. It is already antici- 
pated that a reduction in the data input for successive runs is both 
feasible and desirable. 

Although fourth order Runge-Kutta numerical integration with auto- 
matic error checking and interval modification is accurate, stable and 
relatively rapid, it is recommended that other integration methods in- 
cluding predictor-corrector types, and third and fifth order Runge- 
Kutta be evaluated. 

The FUNCTION subprograms developed for use with SIDES consist of a 
few of the nonlinearities encountered in the analysis of feedback con- 
trol systems. Consideration has been, and should be further given to 
the development of differentiators, noise generators and function genera- 
tors, 

SUBROUTINE SIDES is intended to provide a simple but moderately 
flexible means for the solution of differential equations. The inte- 
gration technique provides a high degree of accuracy and stability and 
the prepared coding forms simplify the programmer’s task. It is be- 
lieved, therefor, that SUBROUTINE SIDES is a valuable and essential 
supplement to the analog computer. 



14 



9 . Acknowledgement . 

The authors wish to thank the personnel of the Computer Facility, 
U, S, Naval Postgraduate School, who assisted in the development of 
this thesis. Particularly, we wish to thank Dr. J. R. Ward for his 
continued interest and counsel as thesis advisor. 



15 



BIBLIOGRAPHY 



16 



1. Anderson, W. H. The solution of simultaneous ordinary differential 

equations using a general purpose digital computer. Communications 
of the Association of Computing Machinery, June 1960: 355-360. 

2. Hurley, J, R. and Skiles, J, J. DYSAC: A digitally simulated 

analog computer. Proc. of the Spring Joint Computer Conference 
(AIEE) 1963: 69-82. 

3. Milne, W. E. Numerical solution of differential equations. 

John Wiley & Sons., 1953. 

4. Nordsieck, A. On numerical integration of ordinary differential 
equations. Mathematics of Computation, v. 16, no. 77, 

January 1962: 22-49. 

5. Stein, M. L., Rose, J. and Parker, D, B. A compiler with an 
analog-oriented input language, Proc. of the Western Joint 
Computer Conference, San Francisco. 1959. 

6. Stein, M. L. and Rose, J. Changing from analog to digital 

programming by digital techniques. Journal of the Association 
of Computer Machinery. January 1960: 10-23. 

7. Thaler, G. J. and Pastel, M. P. Analysis and design of nonlinear 
feedback control systems. McGraw-Hill Book Co,, Inc, 1962, 



/ 



17 




I. 



APPENDIX A 




I 



1 




SUBROUTINE SIDES PROGRAMMING INSTRUCTIONS 



SUBROUTINE SIDES is intended to facilitate the solution of a class 
of problems involving linear and nonlinear differential equations. It 
was written in FORTRAN 60 for the CDC 1604 Digital Computer. The sub- 
routine will normally be called from the FORTRAN Library tape by a 
program which is being compiled and executed under monitor control. 

This calling program and the associated seven data sections required 
by SIDES for the solution of the problem are supplied on 80 column 
punched cards. Utilization of the SIDES Coding Sheets will greatly 
simplify the programmer's task and serves as a check that problem data 
and control commands are properly positioned and in the correct order. 

A description of the functional sections of this form and instructions 
for their use are given below. In what follows, reference should be 
made to the copy of the coding sheets at the end of this appendix. 

These have been filled out for a typical problem and the resulting 
computer output is also appended. 

1. The Calling Program. 

This section contains the problem equation statements which will 
probably cause the inexperienced programmer the most difficulty. 

The first card contains the job identity, normally just the pro- 
grammer's name and his computer facility box number. 

The second card contains the program name. Any combination of 
seven letters or less is acceptable. 

The next three cards are standard and must be punched as shown on 
the coding sheet. 
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Following these cards is the statement of the differential equa- 
tions to be solved. These equations must be written in the form of n 
simultaneous first order equations, where n must not exceed 30. Thus, 

XDOT(l) = fj(X(l),X(2),...,X(n),T) 

• ••••««. 

XDOT(n) = f^(X(l),X(2),...,X(n),T) 

where, 

T = the independent variable 

X(j) = the jth dependent variable 

XDOT(j) = the first derivative of X(j) with 
respect to T 

X(j) may be identified as the output of the jth integrator in the 
corresponding analog computer setup, while XDOT(j) may be identified as 
the input to that integrator. Notice, however, that the sign reversals 
inherent in the analog operations do not occur here. 

It is often desirable to define other functions of the above vari- 
ables, either so that the functions, f^, can be built up in more than 
one step, or so that a quantity other than one of the variables can be 
printed or graphed. Although any acceptable FORTRAN variable name may 
be used to define these functions, they must be equated to an X(j), 
n< j ^ 48, if they are desired as output. 

Finally, if several solutions are contemplated with changes in some 
parameters, the constants C(l) through C(49) may be introduced. Their 
values may then be entered later on data cards. The constant C(50) 
must never be used except as in the calling program, above. 
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Block diagram manipulations are often useful in the formulation 
of the desired first order differential equations. The following 
example shows how a simple pole-zero transfer function may be reformed 
and expressed in FORTRAN equation statements. 




Note ; 



s+6 







R = C(1)*T 

ERROR = R - X(2) 

X(2) = ERROR - X(l) (permits print-out or graph 

of C) 

XDOT(l) = 3.0*(ERR0R - 2,0*X(1)) 

X(3) = ERROR (permits print-out or graph 

of the ERROR) 

Any legal FORTRAN statement or technique can be used in the call- 
ing program, provided that a transfer of control does not prevent the 
processing of all relevant statements during each step of numerical 
integration. Improper transfer of control has proven to be one of the 
most prevalent types of programming errors made by users of SIDES. In 
the variable step size mode of integration, the equation statements are 
frequently evaluated with too large a step size and consequently too 
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large a truncation error. These values are then automatically reject- 
ed by the program and the equations are re-evaluated with a smaller 
step size. This process is repeated until the truncation error is less 
than the prescribed amount. Care must therefor be taken to ensure that 
values which have been set after a transfer of control on an unaccept- 
able pass are reset correctly on the accepted pass. 

The COMMON statement must not be used unless the programmer is 
completely familiar with the SIDES program. All of the FORTRAN Library 
functions are available for inclusion in the equation statements. Addi 
tional FUNCTION subprograms applicable to common control system non- 
linearities have been developed for use in conjunction with SIDES 
and will be discussed later. If the user desires to write his own 
FUNCTIONS or SUBROUTINES, they may be added to the calling program in 
the conventional manner. Comments are permissible throughout the call- 
ing program if a C is placed in column one. 

The calling program must be terminated by: 

GO TO 1 

END 

END 

blank card 

2. The Data and Option Cards. 

C Values . One of the options in parenthesis must be lined out. 

The values of the coefficients, C(j), not simply zeroed or held, are 
placed, with decimal point, in any order, on the next three lines. 

Each coefficient is preceded by its subscript, j, right justified, in 
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the appropriate block. Unused lines are lined out and the total 
number of lines (0 — ^3) containing values is noted on the option 
line for every run. 

X Values . All but one of the options must be lined out. On 
any run after the first, either the initial or final conditions of 
the preceding run may be held and automatically set as the initial 
conditions of the new run. The values are entered in the same manner 
as the C values, lining out the unused lines. 

Time . On the initial run, the initial and final problem times 
are entered, with decimal point, and the HELD command is lined out. 

On subsequent runs either option may be chosen. The final time has 
an upper limit of 10,000 time units. 

Integration Step Size . All but one of the options must be lined 

out. 

If the step size(s) are given, their value (s) are entered, with 
decimal points, on the next line with the corresponding termination 
time(s). The last specified termination time must be the same as the 
final time value entered in the preceding section. Line out the next 
two lines. 

If the step size is to be computed, the line for step size data is 

N 

lined out and the order of magnitude, N, (X(j) ^ 10 , -9$N^99) of 
each dependent variable, X(J), l<j<n, is entered, right justified, 
following the appropriate variable subscript on the next two lines. 
Unused variable subscripts should be lined out. If only one dependent 
variable appears in the problem, the second line is lined out. 
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If the HELD command is given, the following three lines are lined 

out . 

Data Print-out . All but one of the options must be lined out. 

If the GIVEN command is selected, the number of integration steps 
between print-outs must be entered in the appropriate block. Although 
a print-out at every 20th step will usually be sufficient, it should 
be kept in mind that the program will automatically stop at an upper 
limit of 300 lines of print-out. The next line must contain the sub- 
scripts, right justified, of the variables to be printed. For print- 
out, the independent variable, T, is specified by the subscript 50 and 
the integration step size by the subscript 49. On the next two lines 
enter the data print-out column headings in the appropriate block under 
the corresponding variable subscript. All three data cards must be 
present when this option is chosen. 

If the HELD or NONE commands are chosen, the following three lines 
are lined out. 

Graph Output . All but one of the options must be lined out. 

If graph output is needed, enter the desired number of graphs or 
curves in the appropriate block. On the next line enter the number of 
integration steps between plotting points. This should be kept small 
to ensure smoothness, remembering, however, that the program will auto- 
matically stop at an upper limit of 900 plotting points per curve. On 
the last line, enter in the designated blocks the variable subscripts 
of the X and Y ordinates, respectively, for graphs (curves) A, B, C, 

D, and E. Unused graph (curve) letters may be lined out. Since the 
first curve on the multi-curve graph sets the scale, the ’’largest** 
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curve must always be specified first. 

If the HELD or NONE commands are chosen^ the last two lines are 
lined out. 

3. Special FUNCTION Subprograms. 

Certain FUNCTION subprograms have been developed for use in con- 
junction with the simulation of linear and nonlinear feedback control 
systems by SUBROUTINE SIDES. These subprograms are called into use by 
the normal FORTRAN techniques and may be used with either integration 
method of SIDES, The arguments may be expressed by variable names or 
numerical values, with decimal point. A brief description of each 
FUNCTION and its arguments follows. 

FUNCTION DELAYl (PELT. VAR. T ) 

FUNCTION DELAYl will provide as output the designated variable's 
value, delayed in time by a prescribed amount. The output of the de- 
lay will be the initial value of the variable until such time as T 
exceeds DELT. Since this FUNCTION involves storage of past informa- 
tion, it can only be used once in a given program. For this reason, 
additional subprograms, FUNCTION DELAY2 and FUNCTION DELAYS are also 
available. The arguments of each of these subprograms are: 

DELT The value of delay time. The length of delay must 

be greater than the integration step size and less 
than the sum of the previous 499 time increments or 
an error print-out will be given and the program 
will stop. 
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VAR The variable name of the quantity being delayed. 

T The independent variable, time. It must be represent 

ed by T in the calling statement. 



FUNCTION DZONE (XINPUT. HALFD. GAIN ) 



FUNCTION DZONE simulates dead zone or free play. 




IN 



GAIN 



Its arguments are: 

XINPUT The variable name of the input quantity. 

HALFD The absolute value of one half of the dead zone 

or band. 

GAIN The value of the gain of the nonlinearity. 

FUNCTION REL2 (XINPUT. OUTPUT ) 



FUNCTION REL2 simulates an ideal two position relay. 



OUT 





OUTPUT 

/ 


TN 
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Its arguments are: 

XINPUT The variable name of the input quantity. 

OUTPUT The absolute value of the output of the relay. 

FUNCTION REL2H1 (XINPUT. HALFD. OUTPUT ’) 

FUNCTION REL2H1 simulates a two position relay with hysteresis. 
Since it involves the storage of past information, it can only be used 
once in a given problem. For this reason, additional subprograms, 
FUNCTION REL2H2 and FUNCTION REL2H3 are also available. 





> 


' OUTPUT 




f - 




i > IN 






> 


< A 




^ 'vi-HALFD 



Its arguments are: 

XINPUT The variable name of the input quantity. 

HALFD The absolute value of one half of the hysteresis 

band. 

OUTPUT The absolute value of the output of the relay. 
FUNCTION REL3 (XINPUT. OUTPim 



FUNCTION REL3 simulates an ideal three position relay. 

OUT 

OUTPUT 






^ IN 
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Its arguments are: 



XINPUT The variable name of the input quantity. 

OUTPUT The absolute value of the output of the relay. 

FUNCTION REL3D (XINPUT. HALFD. OUTPUT ) 

FUNCTION REL3D simulates a three position relay with dead zone. 



OUT 



J 


OUTPUT 






^ HALFD 





Its arguments are: 

XINPUT The variable name of the input quantity. 

HALFD The absolute value of one half of the dead zone. 

OUTPUT The absolute value of the output of the relay. 



FUNCTION REL3DH1 (XINPUT. DROPOUT. PULLIN. OUTPUT ) 



FUNCTION REL3DH1 simulates a three position relay with dead zone 

and hysteresis. Since it involves the storage of past Information, it 

can only be used once in a given problem. For this reason, additional 

subprograms, FUNCTION REL3DH2 and FUNCTION REL3DH3 are also available. 

OUT 



OUTPUT 



r' 





' 












\ ^PULLIN 

^ DROPOUT 
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Its arguments are: 

XINPUT The variable name of the input quantity, 

DROPOUT The absolute value of input at which the relay 

will drop out. 

PULLIN The absolute value of input at which the relay 
will pull in. 

OUTPUT The absolute value of the output of the relay. 
FUNCTION SATAMP (XINPUT, CUTOFF, GAIN ) 



FUNCTION SATAMP simulates a saturating amplifier. 




Its arguments are: 

XINPUT The variable name of the input quantity. 

CUTOFF The absolute value of the output when it has 
saturated. 

GAIN The value of the linear gain of the amplifier. 



4. Sample Problem. 

The demonstration problem presented here illustrates the method 
of data presentation required by SIDES. The completed SIDES Coding 
Forms indicate the prescribed method of writing equation statements 
and selecting program options. The actual computer output appears next. 
Both data print-out and single and multiple curve graphs are included in 
this output. 
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••JOB* W.T. DOOR BOX 57 
PROGRAM EXAMPLE 

DIMENSION X(50),XD0T130) ,C(50) 

C(50J=1.0 

1 CALL SIDES (T,X,XDOT,C) 

ERROR-C( 1 ) ~X ( 1 ) 

C FCBCK IS INCLUDED tO DEMONSTRATE USE OF LIBRARY FUNCTIONS 

FCBCK=ABSF (S INF ( xm )) /I 0000^0 
RINPUT=ERROR 

0R0P0UT=C(2) ' 

PULLIN=C(3J 

RELAY=REL30H2 ( R I NPUT,C( 2 ) t PULL IN,C ( 4 )) 

SIGNAL=RELAY+FDBCK ' ' 

XD0Tm=X(2) 

XD0T(2)=9^0*SIGNAL-4.0*XI2) 

X(3)=ERR0R 

X(4I=RELAY 

GO TO 1 

END 

END 
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SAMPLE PROBLEM W.T. DOOR BOX 57 

OUTPUT DATA FOR RUN NUMBER 1 



TIME 



OUTPUT 

POSITION 



.OOOOOE+00 
. 14270E+00 
.34270E+00 
.54270E+00 
.59006E+00 
.65054E+00 
.05O54E+OO 
.10505E+01 
. 12505E+01 
. 14505E+01 



.15477E+01 
. 15490E+01 
. 17282E+C1 
. 17733E+01 
.18238E+01 
.20238E+01 
.22238E+01 
.24238E+01 
.26238E+01 
.28238E+01 



.30238E+01 

.32238E+01 

.34238E+01 

.36238E+01 

•38238E+01 

.40238E+01 

.42238E+01 

•44238E+01 

.46238E+01 

.48238E+01 



OUTPUT 

VELOCITY 



•.500COE + 00 
.37448E+00 
.45759E-01 
.60629E400 
.750COE+00 
.91472E+00 
.12441E+01 
. 13950E+01 
.14620E+01 
. 14921E + 01 



.150C0E+01 
. 150C1E + 01 
.13333E+01 
. 12500E+0I 
.n595E + 0l 
.93691E+00 
.83690E400 
.79198E+00 
.77181E+00 
.76277E+00 



.75872E+00 

•75692E+00 

.75613E+00 

.75579E+00 

.75566E+00 

.75561E+00 

•75561E+00 

.75563E400 

.75565E+00 

.75568E+00 



NORMAL STOP AT FINAL TIME 



GRAPH TITLED •• SAMPLE PROBLEM 



GRAPH TITLED •• SAMPLE PROBLEM 



ERROR 



.IOOOOE+00 
.15244E+01 
.25435E+01 
.30014E401 
.30659E+01 
..2407IE+01 
.1O017E+O1 
•4861 5E+00 
•21857E+00 
•98331E-01 



•66736E-01 
• 48739E-01 
.17029E+01 
■• 19791E401 
.161 72E+01 
••72653E400 
•32636E+00 
•14655E+CO 
••65763E-01 
•29463E-01 



-.131 53E-01 
-.5825CE-C2 
-.25323E-02 
-.10529E-C2 
-.3881 IE-03 
-.89423E-04 
.44782E-04 
.10509E-03 
.13218E-03 
.14436E-03 



RELAY 

OUTPUT 



. 15000E+01 
.13745E+01 
.95424E+00 
.39371E+00 
.25000E+00 
.85283E-01 
.24611E+00 
.39503E400 
•.46198E400 
•.49208E+00 



. 15000E + 01 
. 15000E+C1 
.1500CE+01 
.15000E+01 
.OOOOOE+00 
.OOOOOE+CC 
.OOOOOE+00 
.COOOOE+OO 
.OOOOOE+00 
.OOOOOE+CO 



■.50000E + 00 
.50008E+00 
.33330E400 
.25000E+00 
•. 15953E + 00 
.63092E-01 
. 16310E+00 
.20802E+00 
.22819E+00 
.23723E+00 



.OOOOOE+00 
-. 15000E + C1 
-. 15000E+01 
, .COOOOE+00 
.OOOOOE+00 
.OOOOOE+CO 
.OOOOOE+00 
.OOOOOE+00 
.OOOOOE+OO 
.OOOOOE+OO 



.24128E+00 

.24308E+00 

.24387E+-00 

.24421E+-00 

.24434E+00 

.24439E+00 

.24439E+00 

.24437E+00 

.24435E+00 

.24432E+00 



.OOOOOE+OO 

.OOOOCE+00 

.OOOOOE+OO 

.COOOOE+-00 

.OOOOOE+OO 

.OOOOOE+OO 

.OOOOOE+OO 

.OOOOOE+OO 

.OOOOOE+CO 

.OOOOOE+OO 



W.T. DOOR? BOX 57 RUN 1 



W.T. DOOR BOX 57 RUN 1 



Printout for second run not included. 
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SAMPLE PROBLEM 
2 RUNS ARE CALLED FOR. 



W. T. DOOR BOX 57 



INPUT DATA FOR RUN NUMBER 1 



THE NON-ZERO DATA COEFFICIENTS ARE AS FOLLOWSt 



CC 1) = 
C( 2) » 
C( 3) * 
Cl 4) = 



.lCOE+01 

.250E+00 

.500E+00 

•150E+01 



THE NON-ZERO INITIAL CONDITIONS ARE AS FOLLOWSt 

XOI 1 ) = -.500E + 00 
XOI 2) = .lOOE+00 



THE TIMING DATA ARE AS FOLLOWS* 



INITIAL TIME 
FINAL TIME 



.OOOOOE+00 

.50C00E+01 



THE STEP SIZE IS COMPUTED, BASED ON 
THE ESTIMATED ORDER OF MAGNITUDE OF 
THE VARIABLES AS LISTED BELOW. 



XI 1) = l.OE Q 
XI 2) = l.OE-3 



PRINT SUMMARY 

I ' 

20 INCREMENTS BETWEEN PRINTOUTS. 
THE VARIABLES PRINTED ARE, , 



XI50) TIME 

XI 1) OUTPUT 

XI 2) OUTPUT 

XI 3) ERROR 

XI 4) RELAY 



GRAPH SUMMARY 

2 SINGLE GRAPHS. 

PLOT EVERY POINT. 

GRAPH .A IS XI50) 
• GRAPH B IS XI50) 






POSITION 

VELOCITY, 

OUTPUT 



VS. 

VS. 



XI 

XI 



1) 

3) 



; ' M. ■ 



"t , » ; 






* XI50) REPRESENTS TIME. 
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SAMPLE PROBLEM W.T. DOOR BOX 57 

OUTPUT DATA FOR RUN NUMBER 1 



TIME 



OUTPUT 

POSITION 



OUTPUT 

VELOCITY 



ERROR 



RELAY 

OUTPUT 






.OOOOOE+00 
. 14270E + 00 
.342706+00 
.54270E+00 
.59006E+00 
.65054E+00 
.85054E+00 
.105056+01 
.125056+01 
. 14505E+01 



•.500C0E + 00 
-.374486+00 
.45759E-01 
.60629E+00 
.750C0E+00 
.914726+00 
.124<1E+01 
.139506+01 
.146206+01 
. 14921E+01 



.lOOOOE+00 
.15244E+01 
.254356+01 
.30014E+01 
.306596+01 
.240716+01 
.108 176+01 
.486156+00 
.218576+00 
.983316-01 



.150006+01 
.137456+01 
.954246+00 
.393716+00 
.250006+00 
.852836-01 
-.2461 16+00 
-.395036+00 
-.461986 + 00 
-.492086+00 



. 150006+01. ' 
.150006+01 
.150006+01 
.150006+01 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 



.154776+01 
. 154906+01 
. 172826 + 01 
.177336+01 
.182386+01 
.202386+01 
.222386+01 
.242386+01 
.262386+01 
.282386+01 



.150006+01 
. 150016 + 01 
.133336+01 
.125006+01 
.115956+01 
.936916+00 
.836906+00 
.791986+00 
.771816+00 
.762776+00 



.667366-01 
.487396-01 
-.170296+01 
-.197916+01 
-.161 726 + 01 
-.726536+00 
-.326366+00 
-.146556+00 
-.657636-01 
-.294636-01 



-.500006+00 

-.500086+00 

^.333306+00 

-.250006+00 

-.159536+00 

.630926-01 

.163106+00 

.208026+00 

.228196+00 

.237236+00 



.000006+00 
-. 150006 + 01 
-.150006+01 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 



.302386+0 

.322386+0 

.342386+0 

.362386+0 

.382386+0 

.402386+0 

.422386+0 

.442386+0 

.462386+0 

.482386+0 



.758726+00 

.756926+00. 

.756136+00' 

.755796+00 

.755666+00 

.755616+00 

.755616+00 

.755636+00 

.755656+00 

.755686+00 



.131 536-01 
.582506-02 
.253236-02 
-.105296-02 
-.3881 16-03 
.894 236^04 
.447826-04 
.105096-03 
.132186-03 
.144366-03 



.241286+00 

.243086+00 

.243876+00 

.244216+00 

.244346+00 

.244396+00 

.'244396+00 

.244376+00 

.244356+00 

.244326+00 



.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
.000006+00 
^000006+ 00 



NORMAL STOP AT FINAL TIM6 



GRAPH TITLED .. SAMPLE PROBLEM 



W.T. DOOR BOX 57 RUN 1 



GRAPH TITLED 



SAMPLE PROBLEM 



, W.T.lOOOR BOX 57 RUN 1 



•' I 



Printout for second run not Included, 
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\ • 



GRAPH A 

I 

GRAPH B 
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X AXIS SCALE = ■ 1.00E+00 
Y AXIS SCALE = 5„00E-01 

SAMPLE PROBLEM W.T. DOOR BOX 57 RUN 2 GRAPH M 
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LISTING OF SUBROLTINC SIOFS 



SUBRGUTIOt SIOE.S (T, X, XCCf, C) 

DIMENSION XI 50). XD0T{.^0 3, C J 50 ) , IT;rLUlO», JIITLl ilft), XRC30), 

1 IP(IO), IGi20). X0i30), PRHO), G'U 1 0 ) , XlIOOOJ, CTU), 

2 YU90C), X2(900), Y2t90C5, X3I900), Y3I9CC), XUiOCC), 

3 Y4I900), X5(900), Y5I900), A?31), AK I 4 , 30 S , f J 6 ) , 

4 XSI30), XRRI30S, XAS30), XBI30)» KP(30), SCHECM30). 

5 JP ( 30 ) , VI 1500) (wC ) 500 ) ,'’Yt 30) ,D 140 5 

COMMON XI ,X2 ,X3,X4,X5tYl » Y2 , Y3 , vn , y^, , ICHECK, I TI } LE , J II ILE , 

1 AK,F ,A,XS,XR,XRR,XA,XO,KP, JP,XO, IP,PR,V,W,MY,C 

1000 FORMAT (IHI) 

1010 FORMAT (/> 

1020 FORMAT {//) 

1030 FORMAT (///) 

1040 FORMAT MCA8) 

1050 FORMAT ( 12 , AE, A5 , 1 1 , A8) 

1060 FORMAT (36M EOUATIONS/RUNS DATA CARD INCORRECT.) 

1070 FORMAT (35H NUMBER OF EQUAIIO.NS NOT SPECIFIED.) 

1071 FORMAT (41H MAXIMUM NUMBER OF EQUATIONS Al LOWED IS 30.) 

1080 FORMAT (30H NUMBER OF RUNS NOT SPECIFIED.) 

1090 FORMAT (5X,22HCNE RUN IS CALLED FOR./, 

1 48H — ) 

1100 FORMAT {5X,I1,21H RUNS ARE CALLEC FOR./, 

1 48H 

1110 FORMAT I31H MORE THAN NINE PUNS PEQUESTED.) 

1120 FORMAT { A8, A 3 , A6 , 2A4 , 4A3 , A7 , I 1 , A8 ) 

1121 FORMAT (38H CN C VALUE CARD CROSS OUT ONE OPTION.) 

1130 FORMAT (35H C VALUE CARD INCCRPECT OR MISSING.) 

1140 FORMAT (8( I3,F7.2) 1 

1150 FORMAT (27H INPUT DATA FOR fUiN NU."BER ,I1,///,47H THE NCN-7ET0 OAT 
1A COEFFICIENTS ARE AS FOLLOWS,/) 

1160 FORMAT ( 1 OX , 2HC ( , I 2, 4H ) = ,E8.3) 

1170 FORMAT (45H MAXIMUM ALLOWED COETriCIENT SUBSCRIPT IS 49.) 

1180 FORMAT I5X,4hN0NE> 

1190 FORMAT ( A8 , A 3 , A6 , A4 , A7 , A 3 , A7 , 3A8 , A2 , 1 1 , A8 1 

1191 FORMAT (44H CN THE X VALLE CARD SELECT ONE OPTION ONLY.) 

1200 FORMAT (35H X VALUE CARD MISSING UR INCORRECT.) 

12T0 FORMAT (48H MAXIMUM ALLOWED STATE VARIABLE SUBSCRIPT IS 30.) 

1220 FORMAT (48H THE NON-ZERO INITIAL CONDITIONS ARE AS FOLLOWS,/) 

1230 FORMAT (1 OX , 3HXC ( , 1 2 , 4H) = ,E8.3) 

1240 FORMAT ( A8 , A3 , A8 , A4, F7.3 , A8 , A5 , F7. 3 , A8 , A3, A4 , A8 1 

1250 FORMAT (50H CN THE TIME DATA CARD ONE OPTION NOT CROSSED CUT.) 

1260 FORMAT (29H ERROR ON THE TIME DATA CARD.) 

1261 FORMAT ( A8 , A3, A5 , A8 , A7,3 A8, A6 , A4 , A8 ) 

1270 format (51H select ONE OPTION ONLY CN THE STEP SIZE DATA CARD.) 
1280 FORMAT (11,12,12,14(13,12)1 

1290 FORMAT (77H STEP SIZE MODE AND GIVEN DATA ARE NOT CONSISTENT OR TH 
lERE IS A FORMAT ERROR.) 

1300 FORMAT (48H STEP SIZE OPTION: NOT INDICATED OR IS INCORRECT.) 

1310 FORMAT ( 2 ( 1 1 , 2 A6 , F7 . 3, A6 , A7 , F7 . 3 ) ) 

1320 FORMAT (51H FOR GIVEN STEP SIZE OPHON, VALUE OF DT NOT GIVEN.) 

_ ^ VALUE OF IF MUST BE GIVEN) 

1340 FORMAT (32H IHE TIMING DATA ARE AS FOLLOWS,,// 

5X, 15HINITIAL TIME = ,E11.5,/ 

5X,15HFINAL TIME = ,E11.51 



1330 FORMAT (54H FOR G I VEN _ STEP S 1 /C t'P}IC .i., 

1 

2 _ 

1350 FORMAT (5Xi35HTHE STEP SIZE IS rc*-PUTEO, BASED CN,/ 

1 " ’ “ 

2 

1360 FORMAT 



t j A f J J n I n ^ 1 t- r o i c- 1 o 

5X,36FTHE ESTIMATED CROEP. U- MAGNITUDE OF / 
5X,30HTHE VARIABLES AS LISTED BELOW.//) 

( 10X,2HX( ,I2,4H1 = ,4H1.0E,I2) 



1370 FORMAT (5X.15HSTEP SIZE 



= , El 1.5,1 OH FROM T = ,E11.5,8H TC T = 



A • I j ri .j 

1 ,E11.5) 

1380 FORMAT ( A8 , A 3 , 3 A8, A7 , 1 2, 2A8 , Al , A4 , A6 , A4, A4 1 

1390 FORMAT (39H PRINTOUT CARD IS MISSING OR INCORRECT.) 
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CARO, SELECT CNE OPTION ONLY 



1400 FORMAT (42H CN PRINTOUT 
1410 FORMAT (I6,7ll0» 

1420 FORMAT m,8(A8,I2n 

1430 FORMAT ( A8 , A 3 , I 1 , 2A0 , A3, I 1 , 3A8 , A 5 , A4 , A6 , A4 , A2 ) 

1440 FORMAT ( 38H GRAPH DATA CARD MISSING OR INCORRECT.) 

1450 FORMAT I38H CN GRAPH CARO CHOOSE ONE OPTION ONLY.) 

1460 FORMAT ( AS , A3 , 1 2 , A8 ) 

1470 FORMAT ( 8 ( A2 , 1 2 , I 2 , 1 2 , 12 ) ) 

1480 FORMAT (44H MAXIMUM NUMBER OF CURVES OR GRAPHS IS FIVE.) 

1490 FORMAT (4RH SPACING BETWEEN GRAPH POINTS MUST BE SPECIFIED.) 

1500 FORMAT (46H NUMBER OF CURVES OR GRAPHS MUST BE SPECIFIED.) 

1510 FORMAT (52H PLOT EVERY — TH POINT CARD IS MISSING OR INCORRECT.) 

1520 FORMAT (17H PRINT SUMMARY ,/) 

1530 FORMAT 15X,11HN0 PRINTOUT) 

1540 FORMAT (5X,I2,3CH INCREMENTS BETWEEN PRINTOUTS./ 

1 5X, 26HTHE VARIABLES PRINTED ARE,/) 

(5X,18HPRINT EVERY POINT./ 

5X, 26HTHE VARIABLES PRINTED ARE,/) 

(5X,I2,3CH INCREMENTS BETWEEN PRINTOUTS./ 

5X, 24HTHE VARIABLE PRINTED IS,/) 

(5X,18HPRINT EVERY POINT./ 

1 5X, 24HTHE VARIABLE PRINTED IS,/) 

FORMAT ( 10X,2HX( ,I2,1H) , 5X, A8 , 1 X , A8 ) 

FORMAT ?/,5X,24H* Xl50) REPRESENTS TIME.) 

FORMAT {17H GRAPH SUMMARY ,/) 

FORMAT (5X,1CHN0 GRAPHS.) 

FORMAT {6X,I1,15H SINGLE GRAPHS.) 

FORMAT (6X,17H0NE SINGLE GRAPH.) 

FORMAT (6X,4CH0NE MULTIPLE GRAPH LABELED GRAPH M WITH ,11, 16H CUR 
VES PLOTTED.) 

FORMAT (58H0NE MULTIPLE GRAPH LABELED GRAPH M WITH ONE CURVE PLCIT 



1541 FORMAT 
1 

1542 FORMAT 
1 

1543 FORMAT 



1550 

1560 

1570 

1580 

1590 

1591 
1600 

1 

1601 

1 

1610 

1611 

1620 

1630 

1640 

1650 

1660 

1670 

1680 

1690 

1700 

1710 

1730 

1740 

1741 

1750 

1760 

1770 

1780 

1790 

1800 

1810 

1820 

1830 



1 



EO. ) 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 



(5X,I2,27H INCREMENTS BETWEEN POINTS./) 
(6X,17HPL0T EVERY POINT.) 



VS. X(,I2,1H)) 
VS. XI , 12, IH) ) 
VS. XI , 12, IH) ) 
VS. XI , 12, IH) ) 
VS. XI ,12, IH) ) 
VS. XI , 12, IH) ) 
VS. XI , 12, IH) ) 
VS. X(,I2,1H)) 
VS. XI ,I2,1H) ) 
VS. X(,I2,1H)) 



DATA FOR RUN NUMBER ,11//////) 
I1,1X,30H0F OUTPUT DATA FOR RUN 



M0X,13HGRAPH a is X(,I2,8H) 
( lOX, 13HGRAPH B IS X(,I2,8H) 
( lOX, 13HGRAPH C IS X(,I2,8H) 
I lOX, 13HGRAPH D IS X(,I2,8H) 
( lOX, 13HGRAPH E IS X(,I2,8H) 
( lOX, 13HCURVE A IS X(,I2,8H) 
( lOX, 13HCURVE B IS X(,I2,8H) 
( lOX, 13HCURVE C IS X(,I2,8H) 
( lOX, 13HCURVE D IS X(,I2,8H) 
(10X,13HCURVE E IS X(,I2,8H) 
(4X, 10A8 ) 

(5X,27H0UTPUT 
( 1X,5HPAGE , 

//////// ) 

(3X,7(A8,4X), A8) 

I 101 IX, El 1.5) ) 

I//24H STOP AT 



NUMBER ,11,/ 



I//19H 

I//26H 

I//25H 

I//26H 

(//11H 



300 PRINT LINES) 

STOP AT T = 10,000) 

STOP AT 10,000 INCREMENTS) 
STOP AT 900 GRAPH POINTS) 
NORMAL STOP AT FINAL TIME) 
STOP AT XI,I2,10H) = 10,000) 



l//Mr ilUK «i t f > un I - lu , uuu , 

(30H ALL RUNS HAVE BEEN COMPLETED.) 



IND=C(50) 

GO TO (2000, 361C, 3680, 5020, 3523, 3530) , IND 
2000 DO 2030 J=l, 10 
2030 ITITLE(J)=8H 

DO 2040 J=1 , 16 
2040 JTITLE(J)=8H 

C RK I COEFFICIENTS 

CT( 1 )=0.0 
CT(2 )=0.5 
CT(3)=0.5 
CTI4 ) = 1 .0 



40 



c 

c 

2050 

2060 

2070 

2080 

2090 

2091 

2092 

2100 

2110 

2120 

2130 

2140 

C 

2150 

C 

2151 

2152 

2160 

2170 

2180 

2190 

2200 

C 

2220 

2230 

2240 

2250 

2260 

2270 

2280 

2290 

2300 

C 



2060.2050.2060 

2060.2070.2060 



NRC = 0 

PRINT 1000 
READ TITLE CARD 
READ 1040, ( ITITLEI I ), 1=1 ,5) 

READ NUMBER CF RUNS / EQUATIONS CARD 
READ 1050, NE, I CHECK ( 1 ) , ICHECK { 2 ) , NR 
ICHECK(4)=8H ECUATIO 
ICHECK(5)=8H RUNS 
IF I ICHECK(4)-ICHECKM )) 

IF ( ICHECK(5 )-ICHECKI3)) 

PRINT 1060 
STOP 1 
IF INE) 2080,2080,2090 
PRINT 1070 
STOP 2 

IF (NE-30) 2092,2092,2091 
PRINT 1071 
STOP 3 

PRINT 1730, (ITITLEI n.I=l,5) 

IF (NR-1) 2100,2110,2120 

PRINT 1080 

STOP 4 

PRINT 1090 

GO TO 2130 

PRINT 1100, NR 

PRINT 1030 

IF (NR-9) 2150,2150,2140 
PRINT 1110 
STOP 5 

ENTRY POINT FOR ADDITIONAL RUNS 

NRC=NRC+1 

IN0I=1 

NPAGE=1 

MZ=1 

READ C VALUE CARD 

READ 1 120, ( ICHECKII ) , 1=1 , 1 0 ), NC , ICHECK M 1 ) 

ICHECKn5) = 6HZERCED 

ICHECK! 16)=4HHELD 

ICHECK! 17)=ICH£CK! 1 5 ) + IC HECK ! 1 6 1 

ICHECK! 18)=ICHECK!3)+ICHECK!5) 

IF ! ICHECK! 1 7 )- I CHECK ! 18 ) ) 2 1 52, 2 1 51 , 2 1 52 
PRINT 1121 
STOP 6 

IF ! ICHECK!3)-ICHECK! 15) ) 

IF ! ICHECKI5 l-ICHECK! 16) ) 

PRINT 1130 
STOP 7 

DO 2190 J=l,49 
C!J)=0. 

IF INC) 2220,2260,2220 
READ C VALUES 
DO 2250 J=1,NC 

READ 114C ! JP!K) ,XA! K1,K=1,8) 

DO 2250 K=l,8 

IF !JP!K)-49) 2240,2240,2230 
PRINT 1170 
STOP 8 
ITK=JP!K) 

C!1TK)=XA!K) 

PRINT 1150, NRC 
K=0 

DO 2280 J=l,49 

IF !C!J)) 2270,2280,2270 
PRINT 1160,J,CIJ) 

K=K + 1 
CONTINUE 

IF!K) 2300,2290,2300 
PRINT 1180 
PRINT 1020 
READ X VALUE CARO 



ICHECK ! 3) 



2160,2180,2160 

2170,2200,2170 
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m 



READ 1190, (I CHECK! J) , J = 1 , 1 1 ) ,NI , ICHECK(12) 

ICHECK! 13 )=6HZER0£D 
ICHECKC 14)=7HIC HELD 
ICHECKl 15)=7HFC HELD 
I CHECK! 16) = ICHECK! 13) + [CHECK! 14) 

[CHECK! 17 )=[ CHECK! 1 3 ) + [CHECK ! 1 5 ) 

[CHECK! 18 )=I CHECK! 14 )+ [CHECK ! 15 ) 

[CHECK! 19)=[CHECK! 3) +1CHECKI5) 

[CHECKI20)=[CHECK!3)+[CHECK!7) 

[CHECK! 21 )=[CHECKI5) + I CHECK! 7) 

[F ! [CHECK! 161-lCHECK! 19 ) ) 2301,2303,2301 

2301 IF ! [CHECK! 171-1CHECK120 ) ) 2302,2303,2302 

2302 [F ! [CHECK! 181-1CHECKI21 ) ) 2304,2303,2304 

2303 PR[NT 1 191 
STOP 9 

2304 [F ! [CHECK!3)-1CHECK!13) ) 2330,2310,2330 
2310 DO 2320 J=1,NE 

2320 X0!J)=0. 

GO TO 2380 

2330 [F ! [CHECK! 14)-[CHECK!5) ) 2340,2380,2340 
2340 [F ! [CHECK!15)-[CHECK17) ) 2350,2360,2350 
2350 PRINT 1200 
STOP 10 
2360 MZ=2 

DC 2370 J=1,NE 
2370 XO!J)=X!J) 

2380 [F !N1) 2410,2430,2410 
C READ X VALUES 

2410 DO 2420 J=1,N[ 

READ 1140 IJP!K) ,XA! K) ,K=1,8) 

DO 2420 K=1,8 
[TK = JP !K) 

[F![TK-30) 2420,2420,2411 

2411 PRINT 1210 
STOP 12 

2420 XnC [TK )=XA!K) 

2430 PRINT 1220 
K=0 

DO 2450 J=1,NE 

[F!XO(J)) 2440,2450,2440 
2440 PRINT 1220, J, XO!J) 

K=K+-1 

2450 CONTINUE 

[F!K) 2470,2460,2470 
2460 PRINT 1180 
2470 PRINT 1020 
C READ TIME DATA CARD 

READ 1240, ! [ CHECK ! J ), J = 1 ,4), TO, [CHECK! 5) , [CHECK ! 6) , TF , ! [CHECK I 1) , [ 
1=7,10) 

[CHECK! 1 1 )=8hSTART T[ 

[CHECK! 12)=4HHELD 

[CHECK! 13)=1CHECK! 1 1 )+ [CHECK! 12) 

[CHECK! 14)=[CHECKI3)+[CHECK!9) 

[F ! [CHECK!12)-[CHECK! 14) ) 2490,2480,2490 
2480 PRINT 1250 
STOP 13 

2490 [F ! [CHECK! 12)-!CHECK!9) ) 2500,2491,2500 

2491 TF=F!1) 

TC=F!2) 

GO TO 2530 

2500 [F ! [CHECK! 1 1)-[CHECK!3) ) 2520,2521,2520 

2520 PRINT 1260 
STOP 14 

2521 F!1)=TF 
F!2)=T0 

C READ STEP SIZE CARD 

2530 READ 1261, ! [CHECK ![),[= 1 , 1 1 ) 

[CHECK! 12 )=5HG[VEN 
[CHECK! 13)=8hC0MPUTED 
[CHECK! 14)=4HHELD 
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2540 

2550 

2560 

2570 

C 

2580 

2590 

2600 

C 

2601 

2610 

2620 

2630 

C 

2640 

2650 

2660 

2670 

2680 

2690 

2700 

2710 

2720 

2730 

C 

2740 

C 

2750 

2760 

2770 

C 



I CHECK! 15 1 = 1 CHECK! 12 HIC HECK ( 13) 

I CHECK! 16) = ICHECK! 12 ) + ICHECK! 14) 

I CHECK! 17) = ICHE(;K( 13 ) <- IC HECK ! 14 ) 

I CHECK! 18) = ICHECK! 3) + ICHECK! 6) 

I CHECK! 19)=ICHECK! 3) +ICHECK! 10) 

ICHECK! 20)=ICHECK! 6) +ICHECK! 10) 

IF ! ICHECK! 15)-ICHECK! 18 ) ) 2540,2560,2540 
IF ! ICHECK!16)-ICHECKI 19 ) ) 2550,2560,2550 
IF ! ICHECK! 17)-ICHECK! 20 ) ) 2570,2560,2570 
PRINT 1270 
STOP 15 

IF ! ICHECK! 13 )-ICHECK! 6) ) 2610,2580,2610 

READ COMPUTEC STEP SIZE CARD 

READ 1280, NC, ! NH, KP ! I ) , I = 1 , NE , 2 ) 

IF IND-2) 2590,2600,2590 
PRINT 1290 
STOP 16 

IF INE-1) 2601,2720,2601 

READ COMPUTEC STEP SIZE CARD 

READ 1280, NC, ! NH,KP ! I ) , I=2,NE, 2) 

IF (ND-2) 2590,2720,2590 

IF ! ICHECK! 14)-ICHECK! 10) ) 2620,2720,2620 
IF ! ICHECK! 12 )-ICHECK!3) ) 2630, 2640,2630 
PRINT 1300 
STOP 17 

READ GIVEN STEP SIZE CARD 

READ 1310,ND,NH,NH,F (3), NH,NH,F !4) ,NH,NH,NH, F!5) ,NH,NH,FI 6 ) 
X!49)=F!3) 

IF IND-l) 2650,2660,2650 
PRINT 1290 
STOP 18 

IF !F!3)) 2680,2670,2680 
PRINT 1320 
STOP 19 
NDT=1 

IF IFI4) -TF) 2690,2720,2690 
IF !F!5)) 2700,2670,2700 

NDT = 2 

IF (FI6) -TF) 2710,2720,2710 
PRINT 1330 
STOP 20 

PRINT 1340,TC,TF 
IF IND-2) 2750,2730,2750 
PRINT 1010 
PRINT 1350 

PRINT 1360, IJ,KP! J) , J = 1 ,NE ) 

COMPUTE CRITERIA FOR RK II STEP SIZE 
DO 2740 J=1,NE 

A! J) =10.0**KP(J )«1.0E-04/! TF-TO) 

AI31 )=TF-TO 

INITIAL STEP SIZE SET FOR RK II 

FI3)=A! 31 )»1 .OE-05 

X!49)=F!3) 

GO TO 2770 

PRINT 1370, F!3). TO, F!4) 

GO TO (2770, 2760), NDT 
PRINT 1370, FIS), F!4), F!6) 

PRINT 1020 

READ PRINTOUT CARD 

READ 1380, ( ICHECKII ), 1=1,6), INCPR,!ICHECKIJ),J=7, 13) 

ICHECK! 14)=8HCATA LIS 

ICHECK! 15)=4HHELC 

ICHECK! 16)=4HNCNE 

ICHECK! 17)=ICHECK!3)+ICH£CK! 10) 

ICHECK! 18 )= ICHECK! 3 )+I CHECK! 12) 

ICHECK! 19) = I CHECK! 10) + ICHECK! 12) 

ICHECK! 20 ) = ICHECK! 14 ) + ICHECK! 15) 

ICHECK! 21 )=ICHECK( 14 )+ ICHECK! 16) 

I CHECK! 22) = ICHECK! 1 5 ) ♦ IC HECK ! 1 6 ) 

IF ! ICHECK! 17)-ICHECK!20 ) ) 2780,2800,2780 
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2780 IF ( ICHECKC 18)-ICHECK( 21) ) 2 790,2800,2790 
2790 IF ( ICHECKM9)-ICHECK(22 ) ) 2 810,2800,2810 
2800 PRINT 1400 
STOP 22 

2810 IF i ICHECKIUl-lCHECKl 12 ) ) 2820,2940,2820 
2820 IF ( ICHECKU 5)-ICHECK( 10 )) 2840,2830,2840 
2830 INCPR=IR1 
GO TO 2940 

2840 IF ( ICHECKI 14)-1CHECK( 3) ) 2841,2 050,2841 

2841 PRINT 1390 
STOP 21 

2850 IF (INCPR) 2880,2841,2880 
2880 IR1=INCPR 

C READ X VALUES TC BE PRINTED 

READ 1410, ( IP(K),K=1,8) 

DC 2900 1=1,8 

IF IlPd)) 2900,2890,2900 
2890 NP=I-1 

GO TO 2910 
2900 CONTINUE 

C READ PRINTOUT TITLE CARDS 

2910 READ 1420, NF, ( JT I TLE (J ) , NH , J= 1 , 8 ) 

READ 1420, NH, UT I TLE (K ) , NH , K=9 , 1 6 ) 

C READ GRAPH CARD 

2940 READ 1430, ICHECK (1 ) , ICHECK ( 2 ) , NG , I CHECK ( 3 ) , ICHECK ( 4 ) , ICH ECK ( 5 ) , NH 
1,(ICHECK( I), 1=6,13) 

ICHECK! 14)=8h SINGLE 

ICHECK! 15)=8H MULTIPL 

ICHECK! 16)=4HHELD 

ICHECK! 17)=4HN0NE 

ICHECK! 18)=ICHECK!3)+ICHECK! 6) 

ICHECK! 19)=ICHECK! 3)+ICHECK( 10) 

ICHECK(20)=ICHECK!3)+ICHECK! 12) 

I CHECK! 21 ) = ICHECK! 6) +ICHECK! 10) 

ICHECK! 22)=ICHECK!6)+ICHECK! 12) 

ICHECK! 23)= ICHECK! 10)+ ICHECK! 12) 

ICHECK! 24)= ICHECK! 14 )♦ ICHECK! 15 ) 

ICHECK! 25)=ICHECK! 14 )♦ ICHECK! 16) 

ICHECK! 26)=ICHECK! 14 )+ ICHECK! 17) 

ICHECK! 27)= ICHECK! 15 ) + IC HECK ! 1 6 ) 

ICHECK! 28) = I CHECK! 15 ) + ICHECK! 17) 

ICHECK! 29 )=ICHECK! 16)+ ICHECK! 17) 

IF ! ICHECK! 18I-ICHECKI24 ) ) 2941,2950,2941 

2941 IF ! ICHECK! 19)-ICHECK!25 ) ) 2942,2950,2942 

2942 IF ! ICHECK!2C)-ICHECK!26) ) 2943,2950,2943 

2943 IF ! ICHECK!21 1-ICHECKI27 ) ) 2944,2950,2944 

2944 IF ! ICHECK!22)-ICHECK!28 ) ) 2945,2950,2945 

2945 IF ! ICHECK!23)-ICHECK!29 ) ) 2960,2950,2960 

2950 PRINT 1440 
STOP 24 

2951 PRINT 1440 
STOP 25 

2960 IF ! ICHECK!3 )-ICHECK! 14) ) 2962,2961,2962 

2961 M0D=1 

GO TO 3000 

2962 IF ! ICHECK!6)-ICHECK! 15) ) 2970,2963,2970 

2963 NG=NH 
M0D = 2 

GO TO 3000 

2970 IF !ICHECK(17)-ICHECK!12)) 2980,3130,2980 
2980 IF ! ICHECK! 16)-ICHECK! 10 ) ) 2951,2990,2951 
2990 NG=IR2 

GO TO 3130 

3000 IR2=NG 

C READ PLOT INCREPENT CARD 

READ 1460, ICHECK! 1 ), ICHECK! 2) , INCGR, ICHECK! 3) 

ICHECK! 4)=8hFL0T EVE 

IF ! ICHECK!! )-ICHECK !4) ) 3001,3010,3001 

3001 PRINT 1510 
STOP 26 
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3010 

3011 

3020 

3030 



3040 

3050 

3130 

3140 

3150 

3151 

3152 

3153 

3154 

3155 

3160 

3170 

3180 

3190 

3200 

'3210 

3211 

3212 

3213 

3214 

3215 

3216 

3220 

3230 

3240 

3250 

3260 

3270 

3290 



3300 

3301 

3302 

3303 

3304 

3305 

3306 

3310 



IFJNG) 3020,3011,3020 
PRINT 1500 
STOP 28 

IF(NG-5) 3040,3040,3030 
PRINT 1480 
STOP 26 

READ X VALUES TC BE PLOTTED 
NG2=NG*2 

READ 1470, (NH,IG{K) ,NH, IGIK+1),NH,K=1,NC2,2) 

IF (INCGR) 3050,3050,3130 
PRINT 1490 
STOP 27 
PRINT 1520 

IF (INCPR) 3150,3140,3150 
PRINT 1530 
GO TO 3180 

IF (INCPR-1) 3151,3151,3152 

IF INP-1) 3153,3153,3154 

IF (NP-1) 3155,3155.3160 

PRINT 1543 

GO TO 3170 

PRINT 1541 

GO TC 3170 

PRINT 1542, INCPR 

GO TO 3170 

PRINT 1540, INCPR 

PRINT 1550, {IP(J),JTITLE( J) , JTITLEI J+08) ,J=1,NP) 
PRINT 1020 
PRINT 1570 

IF {NGl 3200,3190,3200 
PRINT 1580 
GO TO 3480 

GO TO (3210, 330G) , MOD 
SINGLE GRAPH SECTION PRINTOUT 
IF (NG-1) 3211,3211,3212 
PRINT 1591 
GO TC 3213 
PRINT 1590. NG 

IF (INCGR-1) 3214,3214,3215 
PRINT 1611 
GO TO 3216 
PRINT 1610, INCGR 
PRINT 1620, IG( 1 ) , IG( 2) 

IF (NG“2) 3260,3220,3220 
PRINT 1630,IG(3),IG(4) 

IF (NG-31 3260,3230,3230 
PRINT 1640,IG(5),IG(6) 

IF (NG-41 3260, 32UO, 3240 
PRINT 1650, IG(7), IG(8) 

IF (NG“5) 3260,3250,3250 
PRINT 1660, IG(9) , IG( 10) 

DO 3270 J=1,NG2 

IF nG(J)-5C) 3270,3290,3270 
CONTINUE 
GO TO 3480 
PRINT 1560 
GO TO 3480 

MULTIPLE GRAPH SECTION PRINTOUT 
IF (NG-1) 33C1, 3301, 3302 
PRINT 1601 
GO TO 3303 
PRINT 1600, NG 

IF (INCGR-1) 3304,3304,3305 

PRINT 1611 

GO TC 3306 

PRINT 1610, INCGR 

PRINT 1670, IG(1),IG(2) 

IF (NG-2) 3260,3310,3310 
PRINT 1680, IG(3),IG(4) 

IF (NG-3) 3260,3320,3320 
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3320 

3330 

3340 

3480 



3500 

3510 



3520 

3521 

3522 



C 

3523 



3524 



C 

3525 

3530 

3540 

3550 

3560 

3570 

3580 



3590 

3600 

3601 

3610 



3620 

3630 

3640 

3650 

3660 

3670 

3671 

C 

3680 

3730 



1700, IG(7),IG(8) 
3260,3340,3340 



5 ) 



UTITLEt J),J = 1 
N‘RC 

3500,3510,3500 
UTITLEC J) ,J = 1 , 16) 



PRINT 1690, IG{5J,IG(6) 

IF (NG“4) 3260,3330,3330 

PRINT 

IF (NG“5) 

PRINT 1710, IG(9),IG(10) 

GO TO 3260 
PRINT 1000 
PRINT 1730, 

PRINT 1740, 

IF (INCPR) 

PRINT 1750, 

PRINT 1010 
T = TO 
DT=Fi3) 

DO 3520 J=1,NE 
X( J) = XC( J) 

GO TO 13521,3524) ,MZ 
DO 3522 K=l,30 
DIK)=0. 

MY(K)=1 
C150)=5. 

RETURN 

SET INDEXES FOR SPECIAL FUNCTIONS 
MY( 1 )=2 
MA=MYI2 ) 

MB=NY( 3) 

MC=RY(4) 

MD=MYI5) 

ME=MY(6 ) 

MF=MYI7) 

MG=NY(8) 

MH=NY(9) 

MI=RY( 10) 

LINES=0 
NCPTS=0 
NUMPTS=0 

ENTRY POINT FROR RK I OR RK 1 1 
C(5C)=6. 

RETURN 

IF (INCPR) 3540,3650,3540 
IF (NOPTS) 3550,3600,3550 
IF { XR0DF(N0PTS,50»INCPR ) ) 

IF IXMODFINOPTS, 10»INCPR ) ) 

IF (XRODFINOPTS, INCPR)) 

PRINT 1000 
NPAGE=NPAGE+1 
PRINT 1741, NPAGE,NRC 
PRINT 1750, (JT ITLEU),J = 1 
PRINT 1010 
GO TO (3601 ,361C),RA 
C(50)=2.0 
RETURN 

LINES=LINES+ 1 
DO 3640 J=1,NP 

IF (IP(J)-5C) 3630,3620,3630 
PR( J)=T 
GO TO 3640 
IPJ=IP( J) 

PR( J)=X( IPJ ) 

CONTINUE 

PRINT 1760, IPRU),J = 1 ,NP) 

IF (NG) 3660,3850,3660 

IF (XNODF(NOPTS, INCGR) ) 3850,3670,3850 

GO TO (3671, 3680), RA 

C(5C)=3.0 

RETURN 

GRAPH POINTS STORAGE 
DC 3750 J=1,NG2 

IF (IG(J)-5C) 3740,3730,3740 
GRU) = T 



3560,3580,3560 

3570,3590,3570 

3650,3600,3650 



16) 
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GC TG 375C 
37»40 IGJ = IGU) 

GRC JJ=X( IGJ ) 

3750 CGNriNUE 
3840 NUMPTS = NUN*PTS+1 
Xl(NUPPTSJ=GRt 1 ) 

YHNUMPrS)=GRC2) 

X2CNUP'PTS )=GR( 3) 

Y2«NUMPTS)=GR(4) 

X3( NUf^PTS)=GR{5) 

Y3« NUf^PTS )=GR{ 6 ) 

X4CNUIXPTS )=GR(7) 

Y4CNUPPTS)=GR(8) 

X5{NUPPTS)=GRC9) 

Y5(NUPPTS J =GR( 1C) 

3850 N0PTS=N0PTS+ 1 

IF (LINES-30C) 3852,3851,3851 

3851 PRINT 1770 
GG TO 8000 

3852 IF (T-loE+04) 3870,3860,3860 
3860 PRINT 1780 

GC TG 8000 

3870 IF (NCPTS-IOCOC) 3890,3880,3880 
3880 PRINT 1790 
GO TO 8000 

3890 IF (KUMPTS-9C0) 3910,3900,3900 
3900 PRINT 1800 
GC TO 8000 

3910 IF IT-IF) 3930,3920,3920 
3920 PRINT 1810 

C EXIT TO GRAPH PLOTTING ROUTINE 

GC TO 8000 
3930 DO 3950 J=l,49 

IF f ADSFUC J) )“l.E+04) 3950,3940,3940 
3940 PRINT 1820, J 

GC TO 8000 
3950 CONTINUE 

C SPECIAL FUNCTION SECTION 

GG TO « 4160, 3960) , PA 
3960 GO TO 14090, 3970), MR 
3970 GO TO C4010,398C),MC 
C DELAY! 

3980 IF iMYIll)“5C0) 4000,3990,3990 
3990 MY{ 1 1 )=0 
4000 MYf 1 1 )=MY( 1 1 )+l 
MX=MY( 1 1 1 
V(MX)=0(3) 

W(MX)=D{ 1 ) 

MY{ 14)=MYI 15) 

MYI 12)=MY{ 13 ) 

4010 GC TO 14050,4020) ,MD 
C DELAY2 

4020 IF (MYn6)-lC00) 4 040,40 30,4030 
4030 MY(16)=500 
4040 MY{ 16)=MYC 16 )+l 
MX=MY{ 16) 

VCMX J=0t7l 
W(MX )=DC5) 

MY< 19)=MYI 20 ) 

MYC 17)=MYI 18) 

4050 GG TO 14090, 4060), ME 
C DELAY3 

4060 IF {MYC21 )“1500) 4080,4070,4070 
4070 M.YI21) = 1000 
4080 MYI21 )=MYi21 )+ 1 
MX=MYC21 ) 

VIMX )=on 1 5 

W(MX )=0I9) 

MYI 24) = MY( 25 ) 

MY122)=MY( 23) 
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nnoooonooooon 



4090 

C 

4100 



4110 

C 

4120 



4130 

C 

4140 



4150 

4160 



4170 

C 



5000 



5010 



5020 

5030 

5040 



6000 

6010 

6020 

6030 

6040 

6050 

6060 

C 



GO TO 14110, 4100, f'F 
REL3CH 1,2 ANC 3 
01141=0' 13) 

DM6 )=0n5) 

Dn8) = Dl 17) 

GO TO 14130,4120 ,/^G 
REL2H 1,2 AKC 3 
Dn9)=D120) 

D(21 )=0(22) 

DC23J=0(24 ) 

GO TO «4150,414O,N'H 

OKLASH 

D{25)=D(26» 

MJ=KY(26) 

GO TO 14160, 4150 ,MJ 
D(27)=D126) 

TAD = T 

DC 4170 L=1,NE 
XA{L)=xm 
XBl L)=X(U 

EXIT TO RK I OR RK I I 
GO TO 15000,7000 ,ND 



BEGINING OF RUNGA KUTTA I 



DO 5030 1=1,4 

T=TAD + CTn )»0T 
DC 5010 J=1,NE 

X( J)=XA?J)+CT( I)*AK( I-l , J) 

C(50)=4.C 

RETURK 

DC 5030 J=1,NE 

AK{ I I J )=OT»XDOTU) 

DO 5040 J=l, NE 

XA(J) = XA( J ) -KAKd , J ) +2. »AK( 2, J ) *^2 . « AK ( 3 , J ) ♦■ AK 1 4 , J ) 1 »0. 16 66666667 
GO TO (6000,7010, 7040, 7060, TNOI 

END OF RUNG KUTTA I 

IF (NCT-l) 6020,6010,6020 
DT=F(3) 

GO TO 6050 

IF (NDT-2) 6040,6030,6040 
IF IT“F(4n 6010,6040,6040 
DT=F(5) 

X(49)=0T 
T=TAD+DT 
DC 6060 J=1,NE 
X( J )=XAIJ) 

LEAVE RK I TC PRINT AND/GR PLOT POINT 
GO TO (3530, 3525 ) ,MA 



BEGINNING OF RUNGA KUTTA II 

THIS SUBROUTINE IS BASED ON THE METHOD DESCRIBED IN COMM. OF 
THE AoC.M., JUNE 1960 (PP. 355 -360). IT CHECKS THE 
TRUNCATION ERROR AT EACH STEP OF THE INTEGRATION AND ADJUSTS 
THE STEP SIZE ACCORDINGLY. THE ACTUAL RUNGE-KUTTA INTEGRATION 
IS PERFC9DE4 BY RKUTTA, WHICH IS CONTAINED IN THIS 
PROGRAM (5000 SERIES). THE ARGUMENTS ARE, 

NE= NUMBER OF (FIRST ORDER) EQUATIONS. MAXIMUM NE=30 
T = TIME AT START OF INTEGRATION STEP (UPDATED BY THIS 
SUBROUTINE AFTER THE COMPLETION OF EACH STEP). 

DT = STEP SIZE AT START OF EACH STEP (ALSC UPDATED). 

Xm = THE N DEPENDENT VARIABLES (ALSO UPDATED). 

A(I) = THE SPECIFIED ALLOWABLE ERROR PER UNIT TIME FOR EACH 
OF THE DEPENDENT VARIABLES. 
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A? 31 ns USf D TO ENTER I HE TOTAL TIf’L. 

NOTE THAT JF NECESSARY DT IS REDUCED FROM THE VALUE STATED 
IN THE ARGUMENT UNTIL THE SPECIFIED ACCURACY HAS BEEN 
ACHIEVEC. 

7000 PR£VOT=nr 
TS = T 
H = CT 

DT=2.C*h 
IND 1=2 
GC TO 5000 
7010 DC 7020 1=1, NE 
XRI n = XA{ ! ) 

7020 XAm=XBm 

7030 DT=H 

INDI=3 
GO TO 5000 
7040 DO 7050 1=1, NE 
7050 XRRm=XAm 

TAD=TS«^DT 
INDI=4 
GC TO 5000 
7060 DO 7070 1=1, NE 
XSI n = XAU ) 

7070 xm=xBn) 

U2=C.C3 

DC 7100 1 = 1 ,NE 

E21 = IXR{ 1)-aS( I) )*0. 06666666667 

E21R = D{MF| APSFIE21 J, ABSF H AR SF I XS ( 11 ) ♦ AB S F I XR U ) )- 1 .9 9« 

1 ABSFUmn*1.0E-09)1 

THIS TENDS TC PREVENT ROUND OFF ERROR FROM TAKING CONTROL 
U= E21R / iM m6.0*H) 

IF (U-U2) 7100,7100,7090 
7090 U2=U 

7100 xsc n=xsn 1-E21 

IF IU2-1o01 7150,7110,71 1C 
7110 IF (H-AI31 1* 1.0E-09) 7120,7120,7130 
7120 DT=AI 31 1»1 .OE-09 

THIS SETS THE MINIMUM STEP SIZE TO l.OE-09 TIMES THE TOTAL TIME 
GC TO 7200 
7130 DC 7140 1=1, NE 
XAC n=XBI I ) 

7140 XRtn=XRRm 

H=0. 5«H 
TAD=TS 
GO TO 7030 

THIS RECYCLES TEE INTEGRATION IF THE TRUNCATION ERROR IS EXCESSIVE 
7150 IF IU2-Cc03n 7160,7170, 7170 
7160 DT=2.G»H 

GC TO 7180 

7170 OT = SCRTFISGRTFIC.5/U2n»H 
7180 IF IDT-AI 31 1*0.001 1 7200,7200,7190 
7190 DT=A(31 )*0o0Cl 
GC TO 7230 

THIS SETS the maximum STEP SIZE TO l.OE-03 TIMES THE TOTAL TIME 
7200 IF (DT-PREVDTl 7230,7210,7230 
7210 IC0UNT=ICCUNT+1 

IF riC0UNT“3) 7230,7220,7220 
7220 IC0UNT=0 

DT=AC31 1*1 .OE-05 

THE ABOVE IS INTENCEC TO CLEAR A STEP SIZE JAM 
7230 T = TS^^2.0eH 
XC491=DT 
DO 7240 1=1, NE 
7240 xm=xs?n 

THIS UPDATES T AND XU). 

DT IS UPDATED BY STATEMENTS 7120, 7160, 7170 AND 7190 
END OF RUNG A KUTTA II 

LEAVE RK II TC PRINT AND/OR PLOT POINT 
GC TO 13530, 3525 } ,MA 
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C GRAPH PLOTTING RCUriNT ( fiCCC SERIES} 



8000 IF (NGI 8010.8320.8010 

8010 GC TC {8020,6030,8040,8050,8060,8070,8090,8100,81 10) ,NRC 

8020 ITITLEC6}=RH RUN 1 
GO TC 8120 

8030 ITITLEI6)^8H RUN 2 
GC TG 8120 

8040 ITITLE{6)=8H RUN 3 
GO TC 8120 

8050 ITITLE{6)=8H RUN 4 
GO TO 8120 

8060 ITITLE(6)=8H RUN 5 
GO TC 8120 

8070 ITITLEI6)=8H RUN 6 
CC TO 8120 

8090 ITITLE?6J=8H RUN 7 
GC TO 8120 

8100 ITITLE{6)=8H RUN 8 
GO TO 8120 

8110 ITITLE{6)=8H RUN 9 

8120 SFX=0. 

SFY=0. 

MINCFFX=0 

MINGFFY=0 

LABELN0=1 1 

MCDE=0 

N2=0 

N1=0 



8121 



8130 

8131 



1 



8132 

8140 

8141 



1 



8142 

8150 

8151 



8152 

8160 

8161 



1 



8162 
8170 
8171 . 

.8180 

'8190 

8200 

8210 



GO TO {8121,6190, ROC 

SINGLE GRAPH SECTION 

M0DCURV=0 

LABEL=4H 

KCL=1 

ITITLEC7)=8H GRAPH A 
CALL CRAPH2 { NUl^PTS, X 1 , Y 1 , 8 , 
KINOFFY.LABELNO 
GO TO (8132,6320.8220) ,KGL 
IF (NG-2) 8180,8140,8140 
ITITLE{7)=8H GRAPH 8 
CALL GRAPH2 (NUPPTS, X2 ,Y2 , 8, 

nincffy.labelnc 

GO TO { 8142,8320,8250) ,KCL 
IF (NG-3) 81 80,8150, 8150 
ititlec7)=8h graph C 
CALL CRAPH2 ( NUN PTS , X3 , Y 3 , 8 , 
MINCFFY.LABELNC 
GC TO (8152,8320,8280) ,KCL 
IF {NG-4) 8180,8160,8160 
ITITLE(7)=8H GRAPH D 
CALL GRAPH2 ( NUK P T S , X4 , Y4 , 8 , 
PINCFFY,LA8ELN0 
GO TO (81 62, 8320,8310) ,KOL 
IF (NG-5) 8180,8170,8170 
ITITLE(75=8H GRAPH E 
CALL GRAPH2 { NURPTS , X5 ,Y5 , 8 , 
PINOFFY,LABELNC 

GC TO 8320 

MULTIPLE GRAPH SECTION 
ITITLE(7)=8H graph M 
LABEL=4H A 

IF (NG-1) 8200,8200,8210 
MCDCURV=0 
K0L=2 

GO TO 8131 
MCDCURV=1 
K0L = 3 



MOOCURV,LAREL,ITI TLE,SFX,SFY, 
, M0CE,N1 ,N2) 



MODCURV, LABEL, IT ITLE, SFX , SFY, 
,M0DE,N1,N2) 



MCDCURV, LABEL, I T I TLE , SFX , SF Y, 
,M0I)E,N1 ,N2) 



MODCURV, LABEL, IT I TLE, SFX, SFY, 
,M0CE,N1 ,N2) 



MODCURV. LABEL, I T I TLE , SF X , SF Y, 
,MOnE,Nl ,N2) 



MINOFFX, 

MINOFFX, 

MINCFFX, 

MINOFFX, 

MINOFFX, 
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GO TO 8131 
8220 MCDCURV=2 
LABEL=4H R 

IF (NG-2) 8230,8230,8240 
8230 M0DCURV=3 
KCL = 2 

GO TO 8141 
8240 K0L=3 

GO TO 8141 
8250 LABEL=4H C 

IF (NG-3) 8260,8260,8270 
8260 M0DCURV=3 
KCL = 2 

GO TO 8151 
8270 K0L=3 

GO TO 8151 
8280 LABEL=4H D 

IF (NG-4) 8290,8290,8300 
8290 MCDCURV=3 
KCL = 2 

GO TO 8161 
8300 KCL=3 

GO TO 8161 
8310 M0DCURV=3 
LABEL=4H E 
GO TO 8171 
8320 PRINT 1000 

IF (NRC-NR) 2150,8330,8330 
8330 PRINT 1830 
STOP 
END 
END 
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AVAILABLE FUNCTION SUBPROGRAMS 



FUNCTION DELAYl ( DEL T , VA R » T ) 

DIMENSION V( 1500) ,W( 1500 ), MY (30) ,0(40) 

COMMON V,W,/MY,D 
MA=MY( 1 ) 

GO ro ( 1 ,2 ) ,MA 

1 D(11=VAR 
0(2 ) =VAR 
0( 3)=T 
MYC2)=2 
MY(3)=2 
MV(4)=2 
MY( 1 1 )=0 
RETURN 

2 0(1)=VAR 
D(3)=T 
MC=FY{ 12) 

GO (0 (3,9),MC 

3 MY(13)=1 

IF (CELT-T) e,4,4 

4 IF (MY(11)-5CC) 7,5,5 

5 PRINT 6 

6 FORMAT (43H PERIOD OF DELAYl EXCEEDS AVAILABLE MEMORY.) 
STOP 31 

7 0ELAY1=D(2) 

RETURN 

8 MY( 13)=2 

9 ME=MY(14) 

D(4)=T-DELT 

IC0UNT=1 

10 DO 12 K=ME,5C0 

IF (V(K)-D{4)) 12,12,11 

11 MF=K 

MY( 15)=K 
GO TO 15 

12 CONTINUE 



ME=1 

ICOUNT=ICnUNT+ 1 

IF ( ICOUNT-3 ) 1C, 13, 13 

13 PRINT 14 

14 FORMAT (63H CELAYl MALFUNCTION. CHECK THAT OELT EXCEEDS INITIAL ST 
lEP SIZE.) 

STOP 32 

15 IF (MF-1) 16,16,17 

16 11=500 
12=1 

GO TO 18 

17 I1=MF-1 
I2=MF 

18 DELAY 1=V«( I 1 ) MW( I2)-W( II ) )*(D(4)-vn 1 ) ) / (V(I2)-V( I 1 ) ) 

END 

END 



FUNCTION DELAY2 (DEBT, VAR, T) 

DIMENSION V( 1500) ,W( 1500 ) , MY ( 30 ) ,D ( 40 ) 
COMMON V,W,MY,D 
MA=MY( 1 ) 

GO TO n ,2) , MA 
1 D(5)=VAR 
D(6)=VAR 
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D{75=T 
FY( 2 )=2 
MYJ 3)=2 
MY{ 5)=2 
MY{ 1 6)=500 
MY« 19S=501 
RETURN 

2 0(5)=VAR 
D(7)=T 
MC=NYJ 17) 

CO TO {3,9),NC 

3 MYM8) = 1 

IF (DELT-T) 8,4,4 

4 IF (RY? 16)-1C00) 7,5,5 

5 PRINT 6 

6 FORMAT (43H PERIOD OF DELAY2 EXCEEDS AVAILABLE MEMORY.) 

STOP 33 

7 DELAY2=D(6) 

RETURN 

8 MVn8) = 2 

9 ME=MYH9) 

D(8)=T-DELT 

ICCUNT=1 

10 DC 12 K=ME,1CC0 

IF «VIK)-D18)) 12,12,11 

11 MF=K 
MY?20)=K 
GO TO 15 

12 CONTINUE 
ME=501 

ICCUNT=IC0UNT+1 
IF (ICCUNT-3) 10,13,13 

13 PRINT 14 

14 FORMAT 163H CELAY2 MALFUNCTION. CHECK THAT DELT EXCEEDS INITIAL ST 
lEP SIZE.) 

STOP 34 

15 IF (MF- 501) 16,16, 17 

16 11=1000 
12=501 
GO TO 18 

17 I1=NF-1 
I2=MF 

18 DELAY2=W( 1 1 ) ♦ ( W ( I 2 )-W U 1 ) ) * I D ( 8 ) -V ( 1 1 ) ) / ( V ( 1 2 ) -V { 1 1 ) > 

END 

END 



FUNCTION DELAY3 (DELT, VAR, T) 

DIMENSION VC 1500) ,WC 1500 ) ,MY ( 30 ) ,D( 40 ) 
COMMON V,W,MY,D 
MA=MY{ 1 ) 

GO TO n,2),MA 

1 nC9)=VAR 
DC 10)=VAR 
DC1 1 )=T 
MY(2)=2 
MYC3)=2 
MYI6)=2 

MYC 21 ) = 1000 
MYC 24 j=1001 
RETURN 

2 D(9)=VAR 
DC 1 1 ) = T 
MC=MYC22) 

GO TO C3,9),MC 

3 MYC23)=1 

IF (DELT-T) 8,4,4 

4 IF (MYC21 )-1500) 7,5,5 

5 PRINT 6 



53 



6 FCR^'AT {4?H PCRIOD OF f)tLAV3 EXCEEDS AVAILABLE MEf-^CRV.f 
STCP 35 

7 nELAY3=C« 10} 

RETURN 

8 MV123}=2 

9 MC = ‘-'Y{24! 

0H2 )=T-DELr 
ICGUNT=1 

10 DC 12 K=KE,15C0 

IF (ViK)-nM2n 12,12,11 

n .'^F^K 

MYS25)=K 
GO TC 15 

12 CCNTINUE 
ME=1001 

ICCUNT=ICnuNT* I 

IF I ICQUNT-3) 10,13, 13 

13 PRINT 14 

14 FORMAT (63H DELAYS MALFUNCTION. CHECK THAT DELT EXCEEDS INITIAL ST 
lEP SIZE.) 

STOP 36 

15 IF (MF-1001) 16,16,17 

16 11=1500 
12=1001 
GO TO 18 

17 I1=MF-1 
I2 = MF 

18 DELAY3 = W{ 1 1 ) M W ( 1 2 )-W II 1 ) ) « ( 0(1 2 ) -V (1 1 ) ) / I V( 1 2 ) -V ( 1 1 ) ) 

END 

END 



FUNCTION OZONE { X I NPUT ,H AL FD , GA I N ) 
IF (ABSFCXINPUn-HALFD) 1,2,2 

1 DZ0NE=0. 

RETURN 

2 OUT=GAIN*iABSF( XINPUT)~HALFD) 
OZCNE=SIGNF( CUT, X INPUT) 

END 

END 



FUNCTION REL2 ( x I NPUT , CU TPUT ) 
R EL2 =S I GNFI OUTPUT, XI NPUT ) 

END 

END 



FUNCTION REL2H1 { X INPUT, HALFD, OUTPUT ) 
DIMENSION MY(30),0(40) 

COMMON MV,D 
MA=MY(1 ) 

GO TC (1,2 ) , MA 

1 MY(2)=2 
MY( 8)=2 
Dn9)=0UTPU.T 

2 IF (0(19)) 5,3,8 

3 PRINT 4 

4 FORMAT (20H REL2H1 MALFUNCTION.) 

STOP 37 

5 IF (XINPUT-HALFC) 6,7,7 

6 REL2H1=D( 19) 

D(20)=D(19) 

RETURN 

7 REL2H1=0UTPUT 
D(20)=OUTPUT 
RETURN 
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8 IF ( XINPUTMi/'LFU 10,10,9 

9 REL2H1=DC 19) 

0120=0(19) 

RETURN 

10 REL2H1=-0Un’Ur 
D(2O=-0UiPU1 
END 
END 



FUNCTION REL2H2 ( X INPUT, HALFO, OUTPUT ) 
DIMENSION MY(30) ,D{UO) 

COMMON MY,0 
MA=MYI 1 ) 

GO TO n,2),MA 

1 MYI2)=2 
MY«8)=2 
0(21 )=OUTPUr 

2 IF (0(21)) 5,3,8 

3 PRINT 4 

4 FORMAT (20H REL2H2 MALFUNCTION.) 

STOP 38 

5 IF (XINPUT“HALFO) 6,7,7 

6 REL2h2=D(21) 

D!22)=n(21 ) 

RETURN 

7 REL2H2=0UTPUT 
0(22 )=GUTPUT 
RETURN 

8 IF (XINPUOHAlFO) 10,10,9 

9 REL2H2=0S21) 

0(22 )=D(21 5 
RETURN 

10 REL2H2=-nuTPUT 
D(22)=-0UTPUT 
END 
END 



FUNCTION REL2H3 ( X I NPUT, HALFO, OUTPUT ) 
DIMENSION MY(30),D(40) 

COMMON MY,D 
MA=MY( 1 » 

GO TO n ,2 ) , MA 

1 MY(2)=2 
MY(8)=2 
D(23)=0UTPUT 

2 IF (0(23)) 5,3,8 

3 PRINT 4 

4 FORMAT (20H REL2H3 MALFUNCTION.) 

STOP 39 

5 IF ( XINPUT-HALFO ) 6,7,7 

6 REL2H3=D(23) 

D(24)=D(23) 

RETURN 

7 REL2H3=GUTPUT 
D(24)=CUTPUT 
RETURN 

8 IF ( XINPUT^^HALFC ) 10,10,9 

9 REL2H3=D(23) 

0(24)=D(23) 

RETURN 

10 REL2H3=-OUTPUT 
0(24 )=“OUTPUT 
END 
END 
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FlJNCriGN REL 3 i X I NPb ! , OU 1 PU U 
IF (XIMPUn 2 , 1,2 

1 REL 3 = C.O 
RETURN 

2 REL 3 =SIGNFCCU! PUT,XINPUr ) 

END 

END 



FUNCnON REL 3 D ( X I NPUT , H ALFU , nUTPUT) 
IF i AP.SFjXi’NPUT J-HALFCJ 1 , 2,2 

1 RCL3C = 0.0 
RETURN 

2 REL3D=S1GNFI'CUTPUT,X INPUT ) 

END 

END 



FUNCTION RELSDHI C X I NPUT , DR GPCU T , PULL I N , OUTPUT 5 
DIMENSION MYC30),DU05 
COMMON MY,D 
MA=MVf 1 ) 

GO TO n ,2 J , ^■A 

1 MY{2J=2 
MYI7}=2 

2 IF I APSFU INPUn-DRCPOUT ) 3 , 4,4 

3 REL3DM=0o 
Dn3)=0. 

RETURN 

4 IF CAPSFI KINPun-PULlINI 6 , 5,5 

5 Dn3J = SIGNF{CUTPUT,XINPun 
REL3DF1 =01 13 I 

RETURN 

6 REL3CH1=DM4) 

Dn 3 )=nn 4 j 

END 

END 



FUNCTION REL3DH2 IX I NPUT , CRGPGUT , PULL I N , OUTPUT } 
DIMENSION MY|30),DI40) 

COMMON MY,D 
MA = MY? 1 ) 

GO TO n,2) ,MA 

1 MY!2)=2 
MYI 7 )=2 

2 IF { APSFIXINPUTl-OROPGUT S 3,4,4 

3 REL3DH2=0. 

01155=0. 

RETURN 

4 IF ( ABSn XENFUT5-PULLIN5 6,5,5 

5 D«15 )=SIGNF( OUTPUT, XI INPUT 5 
REL3DH2 = D{15 5 

RETURN 

6 REL30H2 =01165 
DM5 5=0{165 
END 

END 



FUNCTION REL3DH3 I X I NPUT , DRCPOU T , PUL L I N, OUTPUT 5 
dimension my >30) ,D? 40 5 
COMMON MY,D 
MA=MY(1 5 
GO TO { 1,25,MA 
1 N<YI 2 S =2 
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,yY{7)=2 

2 IF (ABSFC XINFUn-DRCPCUr ) 3,4,4 

3 REL3DH3=0. 

D(17)=0. 

RETURN 

4 IF UPSF( XINFUT )-PULlJN) 6,5,5 

5 D n 7 ) =S I GNF { CUTPUT , X I NPU T ) 
REL3CH3=DM7) 

RETURN 

6 REL3DH3=D( 18 ) 

0(17 )=D[ 181 
END 

END 



FUNCTION SATANP ( X I NPUT , CUTOFF , G A I N ) 
IF (GAIN«ABSFfXINPUT )-CUTGFF) 1,2,2 

1 SATANP=GAIN»XINPUT 
RETURN 

2 SATANP=SIGNF (CUTOFF, XINPUT ) 

END 

END 
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